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We present Molecular Dynamics (MD), Quasi Harmonic Lattice Dynamics (QHLD) and Energy 
Minimization (EM) calculations for the crystal structure of Ne, Ar, Kr and Xe as a function of 
pressure and temperature. New Lennard-Jones (LJ) parameters are obtained for Ne, Kr and Xe 
to reproduce the experimental pressure dependence of the density. We employ a simple method 
which combines results of QHLD and MD calculations to achieve densities in good agreement with 
experiment from K to melting. Melting is discussed in connection with intrinsic instability of the 
solid as given by the QHLD approximation. 
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I. INTRODUCTION 

Lattice Dynamics (LD) is based on the expansion of 
the potential energy of the crystal in powers of the dis- 
placements of the atoms from a reference structured^. 
The harmonic approximation, in which the expansion is 
truncated at the second order, reduces the many-body 
problem to many exactly soluble one-body problems 
and allows a direct |Computation of all themiodynamic 
functions. HarmonicEJ and quasi- harmonicQ~l3 lattice dy- 
namics (HLD and QHLD) are alternative LD strategies, 
which differ in the choice of the reference structure. For 
HLD this is the minimum potential energy (mechanical 
equilibrium) structure, whereas for QHLD it is the min- 
imum free energy (thermodynamic equilibrium) average 
structure at a given temperature T and pressure p. This 
difference is crucial, because the accuracy of the har- 
monic approximation is controlled by the amplitude of 
the displacements from the reference structure. There- 
fore QHLD might remain accurate in regimes where the 
average structure deviates significantly from that of the 
potential minimum, provided that the amplitude of the 
vibrations around the average does not grow too large. 
Conventional HLD is definitely unusable in these situa- 
tions. 

A straightforward way to assess the accuracy of the 
QHLD approximation is to choose a model system and 
to compare the QHLD calculations with the essentially 
"exact" classical mechanics results of Molep.ular Dynam- 
ics (MD) or Monte Carlo (MC) simulationsO. If the same 
model system and interaction potential are employed for 
both QHLD and the classical simulations, all differences 
between the two sets of results must be due to the ap- 
proximations implicit m the two methods. 

Lacks and RutledgeQ have chosen a monoatomic fee 
lattice with Lennard-Jones (LJ) interactions as a model 
to study the behaviour of the QHLD and MC free ener- 
gies as a function of T. They find that the QHLD free 
energy is reasonably accurate up to Tm/2, and rapidly 
deteriorates above that; here Tm is the melting temper- 



ature. These results clearly indicate that QHLD as such 
is not reliable close to Tm- However, since not all ob- 
servables behave in the same wayQ, a general statement 
is not yet possible. In fact it is conceivable, though un- 
likely, that specific quantities might deteriorate so slowly 
with increasing T as to remain usable up to T„i- There 
are a posteriori hints that this behaviour might hold for 
the crystallographic parameters. In fact, in recent QHLD 
calculations for naphthalenaa, benzeneO and argonE3 we 
obtained structural results in acceptable agreement with 
the experiments in the whole range of existence of the 
solids. Since close to the amplitude of the atomic 
vibrations are large and the harmonic approximation is 
not expected to be accurate, these results require further 
investigation. To check whether the good results for the 
crystallographic structure near Tm were a genuine con- 
sequence of the method, we have performed QHLD and 
MD calculations for the structures of the Ne, Ar, Kr and 
Xe crystals as a function of p and 21 

FoUoissiing Lacks and RutledgeEl, we started this 
projectEJ by using Verlet'sM LJ parameters for argon. 
These parameters are known to give good agreement with 
the experimental properties of solid, liquid ai ud g aseous 
argon if many-body interactions are neglectecEjO. With 
this effective two-body potential we obtainecO excellent 
QHLD results for the density p, at least for very low T. 

Given the excellent results for Ar, we expected simi- 
larly good behaviour for the heavier rare gases, Kr and 
Xe. We tried more-thaii^O different LJ models from stan- 
dard compilationa3il3~EII and, much to our amazement, 
none of them gave sensible results. All QHLD calcula- 
tions with these models gave large discrepancies with the 
experimental densities at low T: these are precisely the 
conditions where QHLD was supposed to be accurate. 
Even more puzzling was the fact that the experimental 
density did not follow any regular trend with the atomic 
number once reduced with the literature LJ parameters. 
Since the experimental data for Ne, Ar, Kr and Xe do 
follow a regular trend when reduced with respect to the 
critical point data [pc, Tm, Pc)t^, we have concluded that 
the origin of the discrepancies is the usage of inappropri- 
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ate LJ parameters. 

Perhaps one could have forecast these kinds-ot^ifh- 
culties, since most of the hterature parameterg3lljil3 are 
fitted to data on the gas, rather than the sohd. The 
gas data yield, in principle, the true interaction potential 
between two isolated atoms, whereas for the solid and 
liquid phases one needs effective two body models. Such 
LJ models incorporate an average of the many body in- 
teractions and, therefore, can only be valid in a limited 
density range. We wish to stress that we do not believe 
that the true interaction potential in the rare gases fol- 
low a LJ curve, nor that the literature parameters are 
"wrong" . We aim to map the experimental behaviour 
of the rare gas solids onto properly chosen LJ models. 
Since there are two adjustable potential parameters, it 
is always possible to obtain exact agreement between ex- 
periments and calculations for at least two p,T points, 
and approximate agreement for a limited range of p,T 
values. For our purposes we require agreement at low 
p and T (the regime where QHLD is accurate), and we 
have therefore decided to adjust the LJ models to low p, 
T data. Different choices may be appropriate for other 
purposes. 



II. METHODS 
A. Molecular Dynamics 

Moi is a simulation method in which the classical me- 
chanical equations of motion of the system are integrated 
numerically. The method avoids any assumption of har- 
monicity of the lattice vibrations, at the cost of neglecting 
all quantum effects. We use the MD results as calibration 
data to investigate the convergence of the QHLD results 
towards the classical limit, and to evaluate the effects of 
the vibrational anharmonicity. 

We performed a sequence of MD simulation runs for 
a specified number of particles N ^ pressure p and tem- 
perature T. The simulations employed N=2b(!) particles, 
in a cube with periodic boundary conditioHS, and, us- 
ing Andersen's isothermal-isobaric methodEJ, described 
an fee crystal in contact with a heat bath and subject 
to a hydrostatic pressure. The equations of njiation were 
integrated using the velocity Verlet algorithmtllE3. Since 
numerous Mll^udies of monoatomic LJ systems have 
been reportecHO'Ej, we omit other details of the present 
simulation. 



B. Quasi Harmonic Lattice Dynamics 

In the QHLD method!"! the Gibbs free energy of the 
crystal is computed as the free energy of an ensemble of 
phonons of frequency vc 



G{p,T) = ^^+pV- 



i i ^ 

(1) 

Here $o is the total potential energy of the crystal in its 
average structure (the electronic ground state energy), 
pV is the pressure- volume term, is the zero- 

point energy, and the last term is the entropic contri- 
bution. The structure as a function of p and T is de- 
termined by minimizing G with respect to the unit cell 
voluiaej—Further details of the QHLD method appear in 
Refs.ry. 

For comparison with the MD calculations, where quan- 
tum effects are neglected, it is appronijiate to consider the 
classical limit (i.e. h ^ Q) of eq. (l)cl: 



Gciassic(p, T) = ^a+pV + keT In (-^ 



(2) 



It should be noticed that, although eq. (2) still depends 
on Planck's constant h, its derivatives, and thus the equi- 
librium crystal structure, do not. The difference between 
Gciassic and G equals the zero point energy at K, and, 
as appropriate for a classical limit, tends to disappear at 
very high T {ksT > hiyt). 

If all vibrational effects are neglected (the static lattice 
approximation), both G and Gciassic reduce to the purely 
mechanical part of the Gibbs free energy, Ginech(p) = 
$0 + pVi which is the free energy in the limit of a large 
atomic mass. An energy minimization (EM) calculation 
for Gmech avoids the determination of the phonon fre- 
quencies and is therefore less time consuming than the 
complete calculation. 



HI. CALCULATIONS 

We have used MD and QHLD to compute, as a func- 
tion of p and T, the structure of a monoatomic fee 
lattice with pairwise additive LJ interactions, $(r) 
4e[((T/r)^^ — (a/r)^]. To compare different atomic species 
we use non-dimensional reduced units, e.g. temperature 
T* — T/e, pressure p* = pa^ /e and density p* — pa^ /m; 
here m is the atomic mass. This is the standard-approach 
in LJ calculations for monoatomic systemsBlliHl^. When 
expressed in reduced units, the MD, Gciassic and Gmoch 
results turn out to be independent of the atomic species, 
because the classical mechanics system has no energy, 
length and mass scales other than e, a and m. On the 
other hand, the experiments and the QHLD results with 
the quantum free energy G (eq. 1) depend on Planck's 
constant h, which in reduced units reads h* = h/Vmale., 
and are therefore dependent on the atomic number Zcj. 
Because h* tends to zero for large Z (as m, a and e all 
increase with Z) , the reduced data should monotonically 
converge towards a classical limit in the sequence Ne, Ar, 
Xe, Kr. 
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IV. RESULTS 
A. Argon density as a function of pressure 

In the calculations for solid Ar we have used Verlet' E 
LJ parameters (table ||). We report in fig. 1 the results 
for the argon density p as a function xtfj^ressure p. The 
measured p vs p ed, 4, 77 and 293 KEiEj are compared 
to the corresponding MD and QHLD results. The main 
body of the figure contains the data for moderate pres- 
sures, p < 2 GPa, which we discuss first. At 4 K, where 
the anharmonicity is negligible, QHLD reproduces the 
low p measurements almost perfectly, while MD, which 
neglects the zero-point expansion on the lattice, over- 
estimates the experimental density. At higher T, where 
anharmonicity is large but quantum effects are less im- 
portant, the opposite situation holds and MD is found 
to be more accurate than QHLD. The difference between 
MD and QHLD is only noticeable at very low p and tend 
to vanish for increasing p. 

The good results of the LJ model for p < 2 GPa do 
not extend to higher pressures, as shown by the inset 
of fig. 1 for the results up to 30 GPa. Since in this 
p range MD and QHLD are barelyr-distinguishable and 
both disagree with the experimentsEa, the failure of the 
calculations must be due to deficiencies of the potential, 
rather than to the calculation methpd. After analyzing 
their high p data, Grimsditch et al E3 have stated that for 
very dense solid argon there is no way to reproduce the 
experimental results with a pair potential and that many- 
body potentials have to be incorporated. Even though it 
is not yet ruled out that an effective pair potential with 
a different functiopal form could work, it is clear that the 
Verlet 's LJ modelO is not usable beyond 2 GPa. 

B. Dependence of the density on the atomic species 

As discussed in the introduction, none of the literatiLca 
LJ models for Ne, Kr and Xe that we have testedO'llaTll 
was found appropriate for our purposes. We have there- 
fore developed new LJ models by fitting the QHLD cal- 
culations to the density measurements as a function of p 
at the lowest available T and in a range of moderate pres- 
sures (0-2 GPa). This is the regime where the harmonic 
approximation is valid and, as shown by the results for 
Ar, effective two-body LJ models are usable. The opti- 
mal LJ parameters, which are reported in table |, have 
values in the ranges oi-fh^ typical literature models. The 
Ne, Kr and Xe datalElS used in the fit are compared 
in fig. 2 with the QHLD results with the optimal LJ 
models. The Ar measurements and calculations at 4 K 
already shown in fig. 1 are also reproduced in fig. 2, 
together with the results of a EM calculation for Gmcch 
(static lattice approximation), which in reduced units is 
independent of the atomic species. 



The good agreement between experiments and calcula- 
tions in fig. 2 has no particular significance, as it merely 
indicates that the fits are good. Of greater importance is 
the fact that the experimental data, when reduced with 
this particular set of e and p values, follow a sensible 
trend with the atomic mass m. The reduced density p* 
is lower for Ne, which has the lighter atom and therefore 
the largest zero-point effects, and then increases with m. 
As the classical limit is approached in the Ne, Ar, Kr, Xe 
sequence, the experimental p* tends to a limiting curve, 
which appears to be very close to the static lattice cal- 
culations for infinite m. The Gmech results, which can 
be obtained very economically, can be an useful approx- 
imation if one wants to estimate the effects of pressure 
without all the complications of a complete QHLD cal- 
culation. 

To clarify the role of the quantum effects on the den- 
sity, we have computed p* {p, T) at zero p and T as a 
function of to, considered as a continuously varying pa- 
rameter. The resulting p* is displayed as a function of 
h* = h/^ma'^e in the inset of fig. 2, together with the 
experimental data at the lowest available T and p. These 
data constitute, in effect, a vertical cross section of fig. 2 
along the p = Q line. Since the p*(0, 0) vs h* curve does 
not depend on any physical dimension, all the experi- 
mental data at p « 0, T w 0, if properly reduced, should 
lie close to the curve. This is what happens with the LJ 
parameters of table |l_ biit-Jiot with the other LJ mod- 
els that we have triedl3'tj~t3. This discrepancy is one of 
the problems with the literature parameters mentioned 
in the introduction. 

The p* vs h* curve, which is found to deviate very little 
from a straight line, exhibits the expected trend towards 
a limiting density with increasing atomic number. The 
h* — limit (infinite to) is p* — 1.092, the correct min- 
imum energy density for an fee lattice of LJ particlesliH. 
For increasing h* (decreasing to), the zero-point effects 
become progressively more important and the lattice ex- 
pands. 

C. Temperature dependence of the density 

We report in fig. 3 the densities as a function of T. 
The reduced experimental data for Xe, Kr, Ar and Ne at 
atmospheric pressure are compared to the MD density 
(which in reduced units does not depend on the atomic 
species) and to the QHLD densities. With increasing T, 
quantum effects become less important and the experi- 
mental density difference among the various species de- 
creases. The density differences, however, are still large 
when the melting temperature T^fP is reached. This ob- 
servation indicates that, especially for neon, quantum ef- 
fects are still significant even at melting. 

Since it ignores all quantum effects, MD cannot dis- 
tinguish between the different atomic species and over- 
estimates the density at T = 0, though only slightly for 
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the heavier atoms. For very small T, MD predicts a 
thermal expansion linear in T, whereas experimentally 
the density is initially independent of T . This is a non- 
classical zero-point effect, which is properly reproduced 
by the QHLD calculations. For large T, where the anhar- 
monicity is large but quantum effects decrease, the MD 
results tend to agree with the experiments better than 
QHLD. As expected, the agreement is much better for 
Ar, Kr and Xe than for Ne, which remains a quantum 
particle even close to the melting point. 

QHLD agrees with the experimental densities near to 
T = (for Ne, Kr and Xe, by construction). As T is 
raised, QHLD predicts a thermal expansion larger than 
that experimentally observedp-because of the neglect of 
the vibrational anharmonicityS. The discrepancy contin- 
ues to increase, until we encounter a limiting temperature 
Tm, different for each atomic species, above which the 
QHLD calculations fails because no minimum of G can 
be found. The QHLD curves in figure 3 are truncated at 
r,„, which is found to be very close to the experimental 
melting temperature. Further discussion of this point is 
postponed to the next section. 

The results shown in fig. 3 confirm those of fig. 1 and 
indicate that QHLD and MD are essentially complemen- 
tary methods, since QHLD is accurate at low T, while 
MD is better at high T. We have found that QHLD calcu- 
lations with Gciassic provide an useful connection between 
the two methods, as shown for Ar in fig. 4, where we 
compare the QHLD, MD, and Gciassic densities, pqhld, 
Pmd and Pciassic- At K the Pciassic and pmd results 
coincide, since both structures correspond to the abso- 
lute minimum of the potential energy. As T is raised, 
^classic progressively diverges from pmd, because anhar- 
monicity is increasing, and converges towards pqhld, as 
the quantum effects tend to decrease. The difference be- 
tween pmd and Pciassic is solely due to anharmonicity, 
whereas the difference between /Oqhld and Pciassic is only 
due to quantum effects. Therefore one can estimate the 
anharmonic corrections as pmd — Pciassic, and the quan- 
tum corrections as /Oqhld ~ Pciassic- Since interactions 
between anharmonic and quantum effects are ignored, 
these are first order estimates of the corrections. 

By adding the anharmonic corrections to pqhld, or, 
equivalently the quantum corrections to pmd, one get a 
"corrected" density, pcorr = Pqhld + Pmd - Pciassic (also 
shown in fig. 4). We find that the corrected densities 
agree with the experiments better than either /5qhld or 
Pmd, with deviations around 0.7 % on the average (2 % in 
the worst case) in the whole range of stability of solid Ne, 
Ar, Kr and Xe. The anharmonic corrections are usually 
small, 1.6 % on the average and 3.5 % in the worst case, 
and do not change very much with the atomic species. 
The quantum corrections are largest for Ne at K (11 
%) and, as expected, decrease in the Ar, Kr, Xe sequence 
(3.4, 1.7 and 1.3 %, respectively). By comparison, the 
deviations between experimental and calculated densities 



that this is indeed the typical accuracy of the MD and 
QHLD calculations themselves. 



D. Thermal expansion mechanism and melting 

Since the curvature of the interatomic potentials de- 
creases with increasing distances, lattice expansion leads 
to smaller phonon frequencies and, consequently, to lower 
zero-point and entropic contributions to the free energy 
G (eq. 1). The thermal expansion of the lattice is driven 
by the competition between the potential energy term <I>o, 
which favours structures close to mechanical equilibrium, 
and the entropic term, which favours expanded struc- 
tures with smaller frequencies, and which becomes pro- 
gressively more important as T increases. Fig. 5 shows 
the Gibbs free energy G{p,T), computed as a function 
of the density p, for argon at p = and for several val- 
ues of T. At T « 0, where the only difference between 
$0 and G is the zero-point energy, the equilibrium den- 
sity is slightly lower than that yielding the minimum of 
$0- The equilibrium density decreases even further with 
increasing T, since the entropic effects increase. Even- 
tually a critical temperature Tm is reached above which 
the minimum of G vanishes altogether: the QHLD model 
become unstable for T > Tm- To clarify this behaviour 
we show in fig. 5 the locus of the minima of G, which, by 
definition, represents G vs p along the curve of thermo- 
dynamic equilibrium. We have computed this curve as a 
parametric function of T, G(0,T) vs p{0,T), following T 
from to the disappearance of the minimum, which for 
Ar occurs at T^^ = 0.651 (T„ = 78 K). 

The idea that melting is connected to an intrinsic in- 
stability of -the solid goes back to Herzfeld and ft)^- 
pert MayerEfl and has been discussed many-jtimescJLj. 
As convincingly argued by Ross and Wol£3, the vari- 
ous instability models consider the solid phase only, and 
therefore cannot describe genuine thermodynamic transi- 
tions between solid and liquid phases of equal Gibbs free 
energy. A correct prediction of the p, T melting curve 
requires an accurate calculation of the free [©^ji^gy in the 
two phases. However, the stability criteriaEjtZl are still 
useful: loss of stability, although not necessarily exactljy 
coincident with melting, cannot be too distant from itLi. 

A formal—jcharacterization of the instability Tm is 
obtainedEJ~Ej by examining directly the conditions for 
the existence and stability of a local minimum of G: 
dG/dV — and d^G/dV^ > 0. Here we are considering 
stability with respect to small volume fluctuations away 
from the volume V{p, T) of thermodynamic equilibrium, 
i.e. we consider p, T and V as independently varying 
variablesLd. By expressing G in terms of the Helmotz 
free energy F, G{p,T) = pV + F{V,T), the equifibrium 
and stability conditions reduce to p = — {dF/dV)j. (the 
state equations) and i^d'^ F / dV'^) ^ = (3t/V > 0, where 



for Van Der Waals {i.e. molecH 
found to be in the 1-5 % rangeH' 



crystals are usually 
I. Our results confirm 



Pi 



-V{dp/dV)rp is the isothermal bulk modulus. It 



is the vanishing bulk modulus or, equivalently, the di- 
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vergence of the isothermal compressibility kx = 1/(3t = 
— (SYJSp)rp /V, which is the origin of the instability at 
TmlHil^Ea. This description is entirely consistent with that 
observed in our QHLD computations. 

In these computations, where G is numerically mini- 
mized by varying the molar volume y in a sequence of 
finite steps, we find that for T < Tm the volume con- 
verges to an equilibrium value, whereas for T > Tm the 
calculation fails because no minimum is found, as shown 
by fig. 5. The QHLD instability temperatures Tm are 
compared in table ^ with the experimental melting tem- 
peratures T^P of the rare gases. In fig. 3 the experi- 
mental points end near Tm^, whereas the QHLD curves 
are truncated at Tm- Notwithstanding the shortcomings 
due to the use of a solid-only melting criteriontj, and to 
the neglect of the anharmonic contributions to the free 
energyEI, the computed Tm follow quite closely the experi- 
mental data. As shown by fig. 3, the QHLD calculations 
account for the quantum effects on melting, since they 
reproduce the observed decrease of the reduced T^fP for 
the lighter atoms. This is a purely quantum effect, as, 
according to the classical law of corresponding states, 
the reduced melting temperature of classicalpLJ particles 
should be independent of the atomic specie£3. 

V. CONCLUSIONS 

We have computed the density of solid Ne, Ar, Xe, and 
Kr through QHLD, MD, and EM methods. In this way 
we have controlled the classical (MD), harmonic (QHLD) 
and static lattice (EM) approximations, i.e. estimated 
rather precisely their effects. Anharmonic effects are neg- 
ligible at low temperatures and even close to the melting 
point they account only for a few percent of the density. 
Quantum effects are only important for very light atoms 
(Ne) and tend to decrease with increasing temperatures. 
Motional effects, which are usually large, decrease for 
large atomic masses, low temperatures and high pres- 
sures. 

We find that QHLD can be used to compute the struc- 
tural parameters, with slowly degrading accuracy, for all 
temperatures up to that where the model becomes un- 
stable. This breakdown temperature turns out to be a 
fair estimate of the melting point. The QHLD accuracy 
for the structural parameters is much better than that 
for the free energy. This result is not an artifact but a 
genuine property of the QHLD method, though we do 
not have yet a fully satisfactory explanation for it. 

QHLD and MD are found to be complementary, rather 
than competing, methods. QHLD, where quantum ef- 
fects are accounted for, but vibrational anharmonicity 
is neglected, is the better method at low temperatures. 
MD, which ignores all quantum effects but correctly han- 
dles large amplitude vibrations, is appropriate for solids 
close to the melting point and for fluids. The differences 
between MD and QHLD, in the range where both are 



valid, are small. Since QHLD is much more efficient 
than MD, we think that a QHLD computation should 
be one of the first steps in testing any proposed potential 
model, even for problems where disorder or large ampli- 
tude motion will eventually require MD or MC (Monte 
Carlo) techniques. If one wants to study the effects of 
pressure, rather than temperature, even a Gmoch energy 
minimization (EM) can be a useful starting point. The 
EM method, which neglects both quantum and anhar- 
monic effects, is quite accurate for heavy molecules at 
low temperatures. 

For a given problem it may be important to choose 
the fastest method among those applicable, as the speed 
differences are very large. For Ar on a fast RS/6000 work- 
station we needed 0.2 sec for a EM calculation, 20 sec for 
a QHLD optimization, and 75 sec for a MD simulation 
with 256 atoms and 1000 timesteps, the minimum for 
barely acceptable equilibration and statistics. Though 
these times obviously depend on the specific problem, 
their ratios are expected to be quite typical. 

It should be noticed that for Ne, Kr and Xe we have fit- 
ted the potential models solely to the experimental den- 
sity as a function of p at very low T. Though no T depen- 
dent data has been used in the fits, the models reproduce 
well the density as a function of T, within the accuracy of 
the calculation methods. These results provide some jus- 
tification for the procedure of neglecting thermal effects 
while fitting a potential model. This approach, which is 
usually adopted because it is very convenient, is found to 
lead to acceptable results. 

Our work shows that one can efficiently and accurately 
predict the structure of a solid phase, through the whole 
range of its thermodynamic stability, by using a suitable 
combination of EM, QHLD and MD methods. 
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TABLE I. Data for the rare gases. The table lists LJ po- 
tential parameters e and a, mass rrfc-reduced Planck constant 
h* — h/V ma^e, and experimentaltJ and computed melting 
temperatures and Tm. LJ parameters for Ar are from 

RefJij; those for Ne, Kr and Xe are from our fits. 
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Ne 


38.5 


2.786 


20.183 
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24.553 


22 


Ar 


119.8 


3.405 


39.948 


0.186 


83.806 


78 


Kr 


159.9 


3.639 


83.8 


0.104 


115.763 


105 


Xe 


220.9 


3.962 


131.3 


0.065 


161.391 


145 
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FIG. 1. Argon doiisijUii p (g/cm'^) vs pressure p (GPa). 
Symbols: experimeiitsEjE3 at 4 K (circles), 77 K (triangles) 
and 293 K (squares). Curves: calculations with MD (broken 
lines) and QHLD (solid lines). The data at 293 K Ate also 
reproduced in the inset, for pressures up to 30 GPaEj. The 
MD and QHLD curves are not distinguishable at the scale of 
the inset. 



FIG. 2. Reduced density p* = pa'^/m vs reduced pressure 
p* = pa'^/t for the rare ga^ep at low tempera|tHres {T ~ 0). 



Experiments at 4 K for NeEj, Ar, Kr and XetiJ: symbols as 
indicated in the figure. QHLD calculations with G and Gmcch: 
solid and broken lines, respectively. The inset displays p* vs 
h* = h/Vma^e at p ^p-O: QHLD calculations at T = (solid 
line) and experimentsH on Ne at 3 K, Ar and Kr at 4 K, and 
Xe at 5 K (symbols). 



FIG. 3. Reduced density p* vs reduced temperature 
T* = r/e for rtiie-|rare [gases at atmospheric orpsjim-e. Exper- 
iments on NeHtel, ArS'Ea, KrHH and XeHOEZl: symbols 
as indicated in the figure. QHLD and MD results: solid and 
broken lines, respectively. 



FIG. 4. Reduced density p* vs reduced t«jr»erature T* for 
Ar at atmospheric pressure. ExperimentsE3~Ej: circles; MD 
calculations, "corrected" density, QHLD results with G and 
with Gciassic: lines as indicated in the figure. 



FIG. 5. The solid curves represent the reduced QHLD 
Gibbs free energy G*{p* = 0,T*) — G/e of Ar, computed 
as a function of the reduced density p*, for several values of 
T* (namely T* = 0.1, 0.2, 0.3 .. . 0.7, 0.8), as indicated near 
each curve. The broken line represents the locus of the min- 
ima of G* , as discussed in the text, followed from T* = to 
the disappearance of the minimum at T* — 0.651. 
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